Crystallization of a quasi-two-dimensional granular fluid 



P.M. Reis*, R.A. Ingale and M.D. Shattuck 

Benjamin Levich Institute, 
The City College of the City University of New York 
140th St. and Convent Av., 
New York NY 10031, USA 

We experimentally investigate the crystallization of a uniformly heated quasi- 2D granular fluid 
as a function of filling fraction. Our experimental results for the Lindemann melting criterion, 
the radial distribution function, the bond order parameter and the statistics of topological changes 
at the particle level are the same as those found in simulations of equilibrium hard disks. This 
direct mapping suggests that the study of equilibrium systems can be effectively applied to study 
non-equilibrium steady states like those found in our driven and dissipative granular system. 



Equilibrium statistical mechanics is generally not ap- 
plicable to systems far from equilibrium where both en- 
ergy input and dissipation mechanisms are present, and 
identifying relevant tools for understanding these systems 
poses a serious challenge to the scientific community . 
Granular materials have become a canonical system to 
explore such ideas since they are inherently dissipative 
due to inter-particle frictional contacts and inelastic col- 
lisions. Granular materials also have far reaching practi- 
cal importance in a number of industries, but often ac- 
cumulated ad-hoc knowledge is the only design tool used 
0. The dissipative nature of grains means that any dy- 
namical study requires energy injection, typically involv- 
ing vibration or shear Q. An important feature of this 
class of systems is that the driving and dissipation mech- 
anisms can be made to balance such that a steady state 
is achieved. Recent investigation of such N on- equilibrium 
Steady States have shown that connections with equilib- 
rium statistical mechanics may provide an useful anal- 
ogy. For example, a single particle on a turbulent air 
flow has been shown to exhibit equilibrium-like dynam- 
ics and the nature of the melting phase transition in 
two-dimensional granular system is consistent with the 
KTHNY scenario for melting of equilibrium 2D crystals 



In our study we have developed an experimental sys- 
tem to generate a vibrated quasi-two dimensional granu- 
lar fluid of stainless steel spheres that is uniformly heated 
(i.e., energy injection is spatially homogeneous). In the 
insets (a) and (b) of Fig. we present two such ex- 
amples of typical non- equilibrium steady states for filling 
fractions <f> = 0.60 and <f> = 0.76, respectively. The first 
((/) = 0.60) is a disordered dense fluid; there is a high 
collisional rate and at long times the particles randomly 
diffuse across the cell. The second (</> = 0.76) is crys- 
tallized with each sphere packed into a hexagonal array 
locked by its six neighbors. 



In this Letter we analyze the fluid-to-crystal transition 
as a function of filling fraction. The aim of our study 
is two- fold. Firstly, we make a quantitative character- 
ization of the structural changes in the granular layer 
across this transition using a number of classic mea- 
sures, namely the Lindemann criterion for melting, the 
radial distribution function, and the bond order param- 
eter. Then we apply the novel concept of shape factor, 
recently introduced by Moucka and Nezbeda |6| , to mea- 
sure in detail the topology of the Voronoi cells across the 
crystallization transition. In parallel, we establish a di- 
rect comparison between the behavior of our experimen- 
tal system and that of simulations of equilibrium hard 
disks and test the extent to which the above quantities, 
commonly used in equilibrium systems, can can be used 
to study a non-equilibrium system such as ours. 

Our experimental apparatus is adapted from a geom- 
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FIG. 1: Lindemanmn ratio, 7 m , v.s. filling fraction, 0, for 
a granular layer vibrated at / = 50 Hz and r = 4. The 
dotted horizontal line is located at 7™ = 0.15. Crystallization 
occurs at <j) s = 0.719. Insets (a) and (b) are representative 
experimental frames in the fluid and crystal phases, at cj) — 0.6 
and 4> = 0.76, respectively. 
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etry introduced by Olafsen and Urbach 0. We inject 
energy into a collection of stainless steel spheres (di- 
ameter D=l. 191mm) through sinusoidal vertical vibra- 
tion with frequency / and dimensionless acceleration, 
T = A(2itf 2 /g), where A is the amplitude of vibration 
and g is the gravitational acceleration. The spheres are 
confined in a fixed volume gap set by a horizontal stain- 
less steel annulus (101.6mm inner diameter) and sand- 
wiched between two glass plates. The thickness of this 
annulus is 1.6 D, thus constraining the system to be quasi- 
2D. The top glass plate is optically flat, but the bottom 
plate is roughened by sand-blasting generating random 
structures from 50/ira to 500/im. Upon vibration the 
rough plate homogeneously randomizes the trajectories 
of the particles. We record the dynamics of the system 
using high speed photography at 840Hz and track the 
particle trajectories in a (15 x 15)mm 2 central region. 

The system is horizontal to minimize gravity induced 
effects such as rolling and compaction. We vary the total 
number of particles in the fixed volume cell over a wide 
range: from a single particle to an hexagonally packed 
crystal. We define the filling fraction of the granular 
layer as (j) = N [D/(2R)\ 2 , where N is the total number 
of spheres, with diameter, D, in a cell of radius R = 
50.8mm. We fix the forcing parameters at / = 50Hz 
and T = 4 and systematically vary the filling fraction 
from 0.2 < </> < 0.8. 

To interpret the qualitative change in behavior be- 
tween dense fluid and crystalline phases, as <p is changed, 
we first measure the Lindemann ratio. For a wide range 
of materials, Lindemann found || that a solid melts when 
the vibrational amplitude of its atoms reaches a criti- 
cal magnitude, typically between 10% and 15%, of the 
inter-atomic spacing. The Lindemann ratio in the vicin- 
ity of crystallization is, 7 m = y ((r — (r)) 2 )/L, where r 
is the positional vector of the particles and L is the bond 
length between (Voronoi) nearest neighbors, correspond- 
ing to the average lattice spacing in the crystal phase. In 
Fig. [IJ 7 m is plotted at high values of (j). In the range 
0.652 < (j) < 0.719 a sharp drop in j m is observed and 
above <fi > 0.719, the Lindemann ratio becomes approxi- 
mately constant at 7 C ~ 0.15. There the system freezes 
at (j) s = 0.719 in excellent agreement with the crystal- 
lization point or solidus point, for equilibrium hard disk 
simulations: 0f m = 0.716 0. 

The Lindemann criterion is empirical and contains lit- 
tle information about the structural configuration. For 
this we calculate the radial distribution function g(r), 
which is a standard way of describing the average struc- 
ture of particulate systems 0. In Fig. Efa) we plot 
curves of g(r) for representative (j). For low filling frac- 
tions (e.g., (j) = 0.5) we observe fluid-like behavior, and 
g(r) is peaked at r/D =1, 2 and 3, as is commonly 
seen in hard sphere simulations 0. At higher (e.g., 
(j) = 0.65), g(r) develops an additional shoulder be- 
low the r/D = 2 peak, which at higher densities (e.g., 



(j) = 0.7 and (j) = 0.72) evolves into a distinct peak lo- 
cated at r/D = y/3, signifying hexagonal packing. To 
each g(r) experimental curve in Fig. EJa), we have su- 
perposed a corresponding (dashed) curve from a Monte 
Carlo simulation of equilibrium hard disks recently re- 
ported by Moucka and Nezbeda 6] , for identical values of 
<p. The agreement between the experimental and numer- 
ical curves is remarkable, implying that our experimental 
non-equilibrium granular fluid is adopting structural con- 
figurations identical to those found in equilibrium hard 
disk systems. The only deviations occur near r/D = 1, 
as seen in the inset of Fig. Efa) for (j) = 0.60. This 
discrepancy is due to the out of plane collisions in our 
experiments leading to apparent particle overlap in pro- 
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FIG. 2: (a) Experimental (solid) and numerical (dashed, ex- 
tracted from |3) curves of the radial distribution functions 
for 5 values of <j). The arrow points in the direction of de- 
creasing cj). Inset: Section of g(r) curve for (j) — 0.6. (b) 
Radial distribution function at contact, g(r = D) v.s. fill- 
ing fraction. The dashed line corresponds to the theoretical 
Carnahan- Starling equation. cj>i and cj) s are the liquidus and 
solidus points, respectively. 
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jection, which would not be possible if the system were 
exactly two dimensional. The amount of overlap is con- 
sistent with our layer thickness of 1.6 D. This deviation is 
seen in the plot of g(r = D) (i.e., at contact) which cor- 
responds to the absolute maximum of g(r) and is shown 
in Fig. |2t>). For low filling fractions and up to <p ~ 0.57, 
g(D) follows the theoretical curve of Carnahan-Starling, 
g cs (D) = [16-7(/>]/[16(l-(/>) 2 ], which is usually assumed 
in the kinetic theory equation of state for granular gases 
[ll|, but g(D) is systematically lower than g cs (D) by 
~ %14. For (j) > 0.57 the deviations from g cs (D) in- 
crease up to = 0.652 where there is a discontinuity in 
the curve's slope. For 0.652 < <fi < 0.719 there is a period 
of slower growth of g(D) with <p. This is consistent with 
the scenario of the existence of a fluid phase (<j) < 0.652), 
intermediate/transition phase (0.652 < (j) < 0.719) and 
crystal phase (</> > 0.719). 

In addition to the development of correlations in 
the particle positions, angular correlations also arise 
as <p is increased 0. We measure these using the 
(global) bond-orientational order parameter ^ lobal = 

\l/MJ2^i l l N i E^i ei6 *H where M is the number of 
particles in the observation window, 6ij is the angle be- 
tween the particles i and j and an arbitrary but fixed 
reference axis, and Ni is the number of nearest neighbors 
of particle z, found using the Voronoi construction [l^ . 
In Fig. 03 we plot the dependence of ^ lobal on (j). The 
value of the bond orientational order parameter tends to 
unity in the crystal phase, but ^ lobal < 1 for a disor- 
dered phase. 

As with g(D), three different regions with the same 
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FIG. 3: Semi- logarithmic plot of the bond-orientational order 
parameter, ipe. The first two lines, I and II, are least squares 
fits of the form ip ~ exp[A0] and line III is a linear fit of the 
form ip ~ A<f). The dashed and solid vertical lines are located 
at — 0.652 and <j) s — 0.719, respectively. Inset: Linear 
version of the plot. 



phase boundaries: (\>\ = 0.652 (liquidus point) and <p s — 
0.719 (solidus point) can be identified in Fig. |3| based 
on the slope of ^6 (</>)• The observed behavior is consis- 
tent with the two-step continuous phase transition ob- 
served during equilibrium 2D crystallization where 
the first transition transforms the isotropic fluid phase 
into an hexatic phase with long range orientational or- 
dering but no positional ordering and the second trans- 
forms the hexatic phase into a crystal with both long 
range orientational and positional order. 

Moucka and Nezbeda |6| have recently introduced the 
concept of shape factor which is a sensitive measure to 
further quantify structural changes in the fluid-to-crystal 
transition in 2D. ( is defined at the particle level, by 
employing Voronoi tessellation, as d = Cf /47rSi, where 
Si is the surface area and C{ the perimeter of the Voronoi 
cell of the i th particle. For circles ( = 1 and ( > 1 
for all other shapes (( = 4/tt ~ 1.273 for square, ( = 
7r/5 tan(7r/5) ~ 1.156 for regular pentagons, and ( = 
6/V37T 2 ~ 1.103 for regular hexagons). Therefore, ( is a 
quantifier of the topology of the Voronoi cells associated 
with the individual particles. 

In Fig. Efa) we present a surface plot of the distribu- 
tion of shape factor, P(C,^>), and vertical cross-sections 
of P(C, 4>) for fixed <p are presented in Fig. Bib). We 
superpose numerical (dashed lines) data of Monte Carlo 
of equilibrium hard disks 6], for the same values of 0, 
and find that our experimental results are in excellent 
agreement with the numerical simulations. At low 0, 
P(£) exhibits a broad and flat maximum; the particles 
are randomly distributed and no specific type of cells 
are formed. As </> is increased, P(C) becomes increas- 
ingly localized around the maximum which progressively 
moves towards lower values of £. Eventually, for (j) > 0.65 
the distribution becomes bimodal and a distinct second 
maximum appears. In the vicinity of the crystallization 
point, (j) s = 0.719, the original maximum for high ( values 
disappears while the low ( maximum rises sharply (cen- 
tered at ( ~ 1.1, the value for regular hexagons). Fig. 
E{a) clearly shows the existence of two distinct classes of 
shapes. 

To quantify these classes we follow the classification 
scheme of the Voronoi cells proposed by Moucka and 
Nezbeda. An important point to note is that the lo- 
cation of the minimum of P(C)> where it exists, is only 
marginally dependent on </>, and we set Cmin — 1.159. 
Class A consists of particles with ( < ( min . Class B 
particles have Cmin < ( < Cn, and Class C have ( > ( u 
where ( u = 1.25. The upper bound, ( u , is set such that 
at the filling fraction for which both maxima of P(C) 
have equal heights ((j) ~ 0.65), the number of particles in 
classes A and B are the same. We plot the boundaries of 
cell classes on the surface plot of P(£, <p) in Fig. Ufa). 

In Fig. 0fc) we present the ^-dependence of the frac- 
tion of particles belonging to each of the Classes A, B 
and C. The nature of the previously mentioned special 
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FIG. 4: (a) Surface plot for the probability distribution functions of shape factor, P(£, </>). The value PDF is given by the 
adjacent color bar. The two horizontal dashed lines located at £ = 1.159 and £ = 1.25 are the boundaries of classes, A, B 
and C of the Voronoi cells, as defined in the text, (b) Experimental (solid) and numerical (dashed, extracted from 0) vertical 
cross-sections of the P(f , <p) distribution along 5 values of <j> The arrow points in the direction of decreasing <j). (c) Fraction of 
particles in the A, B and C classes, as defined in the text, as a function of filling fraction. 



filling fraction values of <pi and <p s , that separate the dis- 
ordered liquid, the intermediate/transition phase and the 
crystal phases, becomes clear under this classification. 
4>i = 0.652 is the point at which Class A and Class B 
occurs in the same proportions (the fraction of Class C 
is negligible at this point). <p s = 0.719 is the point for 
which the fraction of Class B has sharply dropped to zero 
and the granular layer consists almost entirely of particle 
whose Voronoi cells are regular hexagons, i.e., crystal- 
lization occurs. It remains to be shown if the interme- 
diate phase between and (j) s is simply a coexistence 
regions as suggested by the lever-like dependence of the 
fraction of Classes A and B, or this is an hexatic phase 
with algebraically decaying orientational order [5(. One 
would need to perform the experiments with a consid- 
erably larger imaging window to have sufficient spacial 
extension to properly test such scalings. 

In conclusion, we have reported detailed experimen- 
tal measures of structural changes during the crystalliza- 
tion transition in a homogeneously heated granular fluid. 
Our results are in excellent quantitative agreement with 
Monte Carlo simulations for the crystallization of equilib- 
rium hard disks. It is surprising that the particles in our 
granular layer adopt equilibrium-like structural configu- 
rations even though the system is both driven and dissi- 
pative, i.e., far from equilibrium. The equilibrium struc- 
tural configurations for hard disks are usually determined 
by an entropy maximization argument Whether we 
are able to explain the observed phase transitions in our 
system with entropic-like arguments similar to those used 
in hard sphere systems is an important question which 
arises from our study and needs further investigation. 
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